##############################################
## Replication Code: Migration and Terrorism
## Figures 6, A9, A10, and Table A1
##############################################
rm(list = ls())


sink("log/log_MigrationReplication.txt") # start log
begin.time<-Sys.time()


## compile functions 
library(bpNet)
#library(mcmcr)
library(abind)
library(ggplot2)
library(ggstance)
library(coda)
library(igraph)
library(coda)
library(stargazer)

set.seed(123456789)

## ------------------------- ##
##setwd('~/Desktop/PSRMreplicationLP/data/EmpData')
source('code/net_plot.R', chdir = TRUE)    ## load code to make figures

##------- Put the data together for analysis
imm <- read.csv("data/EmpData/immigrationData.csv")
load('data/EmpData/W_flow.RData')
rhoZ1 <- read.csv("data/EmpData/terrorZ1.csv")[,-1]
rhoZ2 <- read.csv("data/EmpData/terrorZ2.csv")[,-1]
tt <- as.numeric(rownames(W_flow[[1]]))
imms <- subset(imm, ccode %in% tt)  ## extract observations that do not have missingness in W

 attach(imms)
 ## the outcome and covariates
 Y <- ln_FGTD 
 n <- length(Y)
 X1 <- cbind(geddes1, geddes2,geddes3,geddes4,geddes5, logGNI, logpop, logarea, GINI, Durable, AggSF, ColdWar, ucdp_type2 ,ucdp_type3, log_iMigrantsBA, LDV1)

 ## country id and year id
 Uid <-ccode
 Tid <-year

 ## put together data for analysis
 dataImmigration <- data.frame(Uid, Tid, Y, X1)

 ## Generate Table A1
 tdataImmigration <- dataImmigration
 names(tdataImmigration) <- c("ccode", "year", 
 	                            "Terrorist attacks (ln)", "Personalist regime", "Military regime",
 	                            "Single-party regime", "Monarchy", "Hybrid regime", 
 	                            "GNI pc (ln)", "Population (ln)", "Area (ln)", "Inequality (GINI)", 
 	                            "Durable regime", "Failed state", "Cold War", "Interstate conflict",
 	                            "ucdp_type3", "Domestic conflict", "Lagged DV")


 stargazer(tdataImmigration[,-c(1, 2)])





 ## ------------- gen W matrix ------------ ## 

 W2 <- array(as.numeric(unlist(W_flow)), dim=c(94, 94, 31))

 ## Get dimension of TSCS 
  N <- dim(W2)[1]
  TT <- dim(W2)[3]

  ## correct the W matrix (## The matrix has the diagnal nonzero)
  ## and standardize the matrix
  W3 <- W2
  for (i in 1:TT){
    w <-  W2[ , , i]
    diag(w) <- 0
    ww <- (w/rowSums(w))
    ww[ww=="NaN"] <- 0
    W3[ , , i] <- ww
    }

## combine the network measure and group-level intercept
rhoZ <- cbind(1, rhoZ1)

## We use the random walk process this time for rho_t process
## then no intercept is allowed, but initial values of rho_0 is needed
int <- c(1,  rep(0, 30)) 
rhoZ3 <- cbind(int, rhoZ1)

## number of iterations and burnin
mcmc <- 35000
burnin <- 5000

## Model rho_t-MLST-FE: only with fixed effects 
out1.6 <- bpNet(data = dataImmigration,  
	              W = W3,          
	              index = c("Uid", "Tid"), 
	              Yname = "Y",
	              Xname = c("geddes1", "geddes2","geddes3","geddes4","geddes5", "logGNI", "logpop", "logarea", "GINI", "Durable", "AggSF", "ColdWar", "ucdp_type2" ,"ucdp_type3", "log_iMigrantsBA", "LDV1"), # "Wmember" is taken out because it is put in Z
	              Zname = NULL,   
	              Aname = NULL,   
                      Contextual = c("ucdp_type3"),  
	              Contextual.effect = "both",  
	              lagY = FALSE,                ## already included in X
	              force = "both",             
	              r = 0,
	              flasso = 0,             
	              rhoZ = as.matrix(rhoZ3),     
	              rho_spec = "rw",               
                      niter = mcmc,
                      burn = burnin)

## Model rho_t-MLST-MF: with factors and Bayesian shrinkage
out1.1 <- bpNet(data = dataImmigration,   
	              W = W3,          
	              index = c("Uid", "Tid"), 
	              Yname = "Y",
	              Xname = c("geddes1", "geddes2","geddes3","geddes4","geddes5", "logGNI", "logpop", "logarea", "GINI", "Durable", "AggSF", "ColdWar", "ucdp_type2" ,"ucdp_type3", "log_iMigrantsBA", "LDV1"), 
	              Zname = NULL,  
	              Aname = NULL,  
                      Contextual = c("ucdp_type3"),  
	              Contextual.effect = "both",  
	              lagY = FALSE,                
	              force = "both",              
	              r = 10,
	              flasso = 1,             
	              rhoZ = as.matrix(rhoZ3),     
	              rho_spec = "rw",                
                      niter = mcmc,
                      burn = burnin)


## Model rho-MLS-MF: time-invariant rho
out1.4 <- bpNet(data = dataImmigration,   
	              W = W3,          
	              index = c("Uid", "Tid"), 
	              Yname = "Y",
	              Xname = c("geddes1", "geddes2","geddes3","geddes4","geddes5", "logGNI", "logpop", "logarea", "GINI", "Durable", "AggSF", "ColdWar", "ucdp_type2" ,"ucdp_type3", "log_iMigrantsBA", "LDV1"), 
	              Zname = NULL,  
	              Aname = NULL,   
                      Contextual = c("ucdp_type3"),  
	              Contextual.effect = "both", 
	              lagY = FALSE,                
	              force = "both",           
	              r = 10,
	              flasso = 1,              
	              rhoZ = as.matrix(rhoZ),      
	              rho_spec = "constant",              
                      niter = mcmc,
                burn = burnin)

##--------------- Simulation ends here  ----------
## ----------------Graph-making begins here --------
########################
## Generate Figure 6
########################
### Figure 6-(b)
out <- out1.1
rho <- rep(0, 31)
TT <- 31
modelRho <- out$Rho
## ----------------------- 1 varying rho_t
datap <- cbind.data.frame(1970:2000, apply(modelRho, 1, mean),
		  apply(modelRho, 1, quantile, 0.025),
		  apply(modelRho, 1, quantile, 0.975),
		  c(rho))
names(datap) <- c("time", "mean", "cil", "ciu", "true")
datap$type <- rep("m", TT)
p <- ggplot(datap) + xlab("Time") + ylab("Rho")
p <- p + geom_line(aes(time, mean, colour = type, size = type, linetype = type))
p <- p + geom_line(aes(time, true), size = 0.5, colour = "red", linetype = "dashed")
p <- p + geom_ribbon(aes(x = time, ymin = cil, 
	           ymax = ciu), fill = "#000000", alpha = 0.2)
set.limits <- c("t", "m", "c")
set.labels <- c("0", "Posterior Mean", "95% CI")
set.colors <- c("red", "grey50", "#00000080")
set.size <- c(1, 1, 3)
set.type <- c("dashed", "dashed", "solid")
p <- p + scale_colour_manual(limits = set.limits,
	                             labels = set.labels,
	                             values =set.colors) +
	 scale_size_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.size) +
	 scale_linetype_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.type) +
	 guides(colour = guide_legend(title=NULL, nrow=1),
	                size = guide_legend(title=NULL, nrow=1),
	                linetype = guide_legend(title=NULL, nrow=1))
p <- p + theme(legend.text = element_text(margin = margin(r = 10, unit = "pt"), size = 12),
	               legend.position = "bottom",
	               axis.title = element_text(size=12),
	               axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0)),
	               axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
	               axis.text = element_text(color="black", size=8),
	               axis.text.x = element_text(size = 8),
	               axis.text.y = element_text(size = 8),
	               plot.title = element_text(size = 15,
	                                         hjust = 0.5,
	                                         face="bold",
	                                         margin = margin(10, 0, 10, 0)))
p <- p + scale_y_continuous(limits = c(-0.031, 0.145))+ xlab(" ") + ylab(" ")
a <- 20
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
p <- p+theme(legend.title = element_text(size=16), legend.text = element_text(size=16), legend.key.size = unit(1, 'cm'))
ggsave('plots/MIG/RhoMF.pdf', p, width = 6, height = 6)

beta1 <- as.mcmc(t(out$beta))
rho <- as.mcmc(t(out$Rho))
rhoA <- as.mcmc(t(out$rhoA))


cat("Model rho_t-MLST-MF: mean \n")
print(mean(datap$mean))
cat("\n")

### Figure 6-(a)
out <- out1.6
### rho
rho <- rep(0, 31)
TT <- 31
modelRho <- out$Rho
## ----------------------- 1 varying rho_t
datap <- cbind.data.frame(1970:2000, apply(modelRho, 1, mean),
		                          apply(modelRho, 1, quantile, 0.025),
		                          apply(modelRho, 1, quantile, 0.975),
		                          c(rho))
names(datap) <- c("time", "mean", "cil", "ciu", "true")
datap$type <- rep("m", TT)
p <- ggplot(datap) + xlab("Time") + ylab("Rho")
p <- p + geom_line(aes(time, mean, colour = type, size = type, linetype = type))
p <- p + geom_line(aes(time, true), size = 0.5, colour = "red", linetype = "dashed")
p <- p + geom_ribbon(aes(x = time, ymin = cil, 
	             ymax = ciu), fill = "#000000", alpha = 0.2)
set.limits <- c("t", "m", "c")
set.labels <- c("0", "Posterior Mean", "95% CI")
set.colors <- c("red", "grey50", "#00000080")
set.size <- c(1, 1, 3)
set.type <- c("dashed", "dashed", "solid")
p <- p + scale_colour_manual(limits = set.limits,
	                             labels = set.labels,
	                             values =set.colors) +
	         scale_size_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.size) +
	         scale_linetype_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.type) +
	         guides(colour = guide_legend(title=NULL, nrow=1),
	                size = guide_legend(title=NULL, nrow=1),
	                linetype = guide_legend(title=NULL, nrow=1))

p <- p + theme(legend.text = element_text(margin = margin(r = 10, unit = "pt"), size = 12),
	               legend.position = "bottom",
	               axis.title = element_text(size=12),
	               axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0)),
	               axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
	               axis.text = element_text(color="black", size=8),
	               axis.text.x = element_text(size = 8),
	               axis.text.y = element_text(size = 8),
	               plot.title = element_text(size = 15,
	                                         hjust = 0.5,
	                                         face="bold",
	                                         margin = margin(10, 0, 10, 0)))
p <- p + scale_y_continuous(limits = c(-0.029, 0.145))+ xlab(" ") + ylab(" ")
a <- 20
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
p <- p+theme(legend.title = element_text(size=16), legend.text = element_text(size=16), legend.key.size = unit(1, 'cm'))
ggsave('plots/MIG/RhoFE.pdf', p, width = 6, height = 6)

cat("Model rho_t-MLST-FE: mean \n")
print(mean(datap$mean))
cat("\n")

### Figure 6-(c)
out <- out1.4
rho <- as.mcmc(t(out$rho0))
In <- rep("rho", dim(rho)[1])
rhop <- data.frame(rho, In)
p <- ggplot(data= rhop , aes(x=In, y=rho))+
    geom_boxplot() +  xlab(" ") + ylab(" ")
p <- p + scale_y_continuous(limits =  c(-0.029, 0.145))
a <- 25
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
###
ggsave("plots/MIG/RhoCon.pdf", p, width = 2, height = 6)

cat("Model rho-MLST-MF: mean \n")
print(mean(rho))
cat("\n")

####################################
## Figure A10: factor selection
####################################
p1 <- omegaPlot(out1.1, oxlim = c(-1, 1), plotlength = 2, labelname = "Factor")
ggsave('plots/MIG/factorsMF.pdf', p1, width = 10, height =5)

#####################################
## Generate Figure A11
## the coefficients of three models
#####################################

paraU1 <- as.mcmc(t(out1.1$beta))
paraG1 <- as.mcmc(t(out1.1$rhoA))
paraU2 <- as.mcmc(t(out1.6$beta))
paraG2 <- as.mcmc(t(out1.6$rhoA))
paraU3 <- as.mcmc(t(out$beta))

nameCC <-   c("Intercept", "Personal Regime", "Military Regime","Single-Party Regime","Monarchy","Hybrid Regime","GNI pc (ln)", "Population (ln)", "Area (ln)", "Inequality (GINI)", "Durable Regime", "Failed State", "Cold War", "Interstate Conflict" ,"Domestic Conflict", "Migrant Inflows (ln)", "Lagged DV", "W*Domestic Conflict")
nameC <-   c("Intercept", "Personal Regime", "Military Regime","Single-Party Regime","Monarchy","Hybrid Regime","GNI pc (ln)", "Population (ln)", "Area (ln)", "Inequality (GINI)", "Durable Regime", "Failed State", "Cold War", "Interstate Conflict" ,"Domestic Conflict", "Migrant Inflows (ln)", "Lagged DV","W*Domestic Conflict", "rho0", "DENSITY")

nn <- length(nameC)
nnn <- length(nameCC)
aa <- rbind(summary(paraU1)[[2]], summary(paraG1)[[2]], summary(paraU2)[[2]], summary(paraG2)[[2]], summary(paraU3)[[2]])
model <- c(rep("Multifactor", nn),rep("Fixed Effects",nn), rep("Constant", nnn))
parameter <- c(nameC, nameC, nameCC)
ddd <- data.frame(parameter, aa, model)
names(ddd) <- c("parameter", "low1", "low2", "Mean", "Upper2", "Upper1", "model")

## put those posteriors with similar scales together
nn3 <-  c(7, 10, 11,12 ,16,20);nn2 <- c(8,  14,15,17,19); nn1 <- c(2:6, 9);nn4 <- c(8,  14,15,17)

## Figure A11 upper-left panel
upa <- ddd[ c(nn1,(nn1+20), (nn1+20+20)),]
CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)
CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                       fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels = c("rho-MF", "rho_t-FE", "rho_t-MF" )) + xlab(expression(beta)) + ylab("")+ggtitle("")

a <- 16
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))

ggsave('plots/MIG/beta1.pdf',CIplot1 , width = 7, height = 6)

## Figure A11 upper-right panel
upa <- ddd[c(nn2,(nn2+20), (nn4+20+20)),]
CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)
CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                      fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels = c("rho-MF", "rho_t-FE", "rho_t-MF" )) + xlab(expression(beta)) + ylab("")+ggtitle("")
a <- 16
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))
CIplot1 <-  CIplot1 +  theme(legend.title = element_text(size=16), legend.text = element_text(size=14))
ggsave('plots/MIG/beta2.pdf',CIplot1 , width = 7, height = 6)

## Figure A11 lower-left panel
rows <-  c(nn3,(nn3+20), (nn3+20+20))
upa <- ddd[rows[-18],]
CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)

CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                      fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels = c("rho-MF", "rho_t-FE", "rho_t-MF" )) + xlab(expression(beta)) + ylab("")+ggtitle("")

a <- 16
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))
CIplot1 <-  CIplot1 +  theme(legend.title = element_text(size=16), legend.text = element_text(size=14))
ggsave('plots/MIG/beta3.pdf',CIplot1 , width = 7, height = 6)

## Figure A11 lower-right panel
pp<- c(nn3,(nn3+20), nn2,(nn2+20), (nn2+40), nn1,(nn1+20), (nn1+40), (nn3+40))
ppp <- pp[pp<60]
upa <- ddd[-ppp,]
CIplot1 <- ggplot(upa, aes(colour = factor(model)))
CIplot1 <- CIplot1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2)

CIplot1 <- CIplot1 + geom_linerangeh(aes(y = factor(parameter),
                                         xmin = low2, xmax = Upper2), lwd =1, position = position_dodgev(height = 1/2))

CIplot1 <- CIplot1 + geom_pointrangeh(aes(y = factor(parameter)
                                          , x = Mean, xmin = low1, xmax = Upper1), lwd =.45,  position = position_dodgev(height = 1/2),
                                       fill = "WHITE")
CIplot1 <- CIplot1 + theme_bw()
CIplot1 <- CIplot1 + theme(legend.position="bottom")
CIplot1 <- CIplot1 + labs(colour = "Model")
CIplot1 <- CIplot1+ scale_colour_discrete(labels =c("rho-MF", "rho_t-FE", "rho_t-MF" )) + xlab(expression(beta)) + ylab("")+ggtitle("")
CIplot1

a <- 12
CIplot1 <- CIplot1 +theme(axis.text.x = element_text( hjust = 1, size=a))
CIplot1 <- CIplot1 +theme(axis.text.y = element_text( hjust = 1, size=a))
CIplot1 <-  CIplot1 +  theme(legend.title = element_text(size=12), legend.text = element_text(size=11))
ggsave('plots/MIG/beta4.pdf',CIplot1 , width = 6, height = 3.5)

## The time-average posterior mean
modelRho <- out1.6$Rho; mean(modelRho)
modelRho <- out1.1$Rho; mean(modelRho)
modelRho <- out1.4$rho0; mean(modelRho)

## Figure A9: examples of the migration network in year 1971, 1991, and 2000
#library(igraph)
#exm <- c(2, 22, 31)
#year <- c(1971, 1991, 2000)
#for (i in 1:3){
#    k <- exm[i]
#   g6 <- graph_from_adjacency_matrix(W3[,,k], weighted=TRUE,mode = "directed")
#   E(g6)$width <- E(g6)$weight
#   V(g6)$label <- rownames(W_flow[[1]])
#   V(g6)$size <- 6
#   l <- layout_with_fr(g6)
#   pdf(paste('plots/MIG/Mig',year[i] , ".pdf", sep=""), width=30, height=30)
#   plot(g6, layout=l, edge.arrow.size=0.1, arrow.width=0.1)
#   dev.off()
#  }


print(Sys.time()-begin.time) 
sink() # close log   

